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Abstract. We present a nonrelativistic calculation of the rotation- vibration levels 
of the molecular ions H^, Dj" and HD + , relying on the diagonalization of the exact 
three-body Hamiltonian in a variational basis. The J = 2 levels are obtained with a 
very high accuracy of 10 -14 a.u. (for most levels) representing an improvement by five 
orders of magnitude over previous calculations. The accuracy is also improved for the 
J = 1 levels of and with respect to earlier works. Moreover, we have computed 
the sensitivities of the energy levels with respect to the mass ratios, allowing these 
levels to be used for metrological purposes. 

cr. 

PACS numbers: 33.15.Pf 

■ 1. Introduction 

In recent years, precise calculations in the hydrogen molecular ion Hg" and its isotopes 
have attracted interest, because these systems appear as promising candidates for the 
metrology of the electron to proton mass ratio Mp/m e , or the ratios of the nuclear 
masses Mp/Mp [H [21 [3]. Almost all the rotation- vibration levels of HD + , H^" and 
including relativistic and radiative corrections have been computed by R. E. Moss 
[HOE] with an accuracy of 10~ 9 a.u. Recent progress in variational calculations has 
allowed to improve the accuracy of the nonrelativistic calculations up to 10~ 14 or even 
10~ 18 for the lowest levels [H [21 El E], while the accuracy on relativistic and QED 
corrections is also improving and should reach 10~ 10 — 10 -11 a.u. [T0l[lT]. However, the 
optical transitions that may be used in metrology experiments also involve more excited 
states [H [2] especially in the perspective of comparing several transition frequencies to 
test the time independence of the mass ratios as proposed in [3]. It would then be useful 
to extend the high-precision calculations as far as possible into the rotation-vibration 
spectrum. 
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To our knowledge, our method relying on the use of perimetric coordinates is the 
only one that allows to maintain a high accuracy up to the most excited vibrational 
states with reasonable numerical means [H [2]. Indeed, this set of coordinates takes 
advantage of the dynamical symmetries of the three body Coulomb problem, so that 
it is possible to choose a basis in which the Hamiltonian has strict coupling rules, and 
whose wavefunctions have the correct long-range behaviour. However, as was already 
discussed in [1], this method becomes less and less advantageous when the value of J 
increases ; only the J = and J = 1 states have been computed so far, with an accuracy 
of 10 -14 a.u. for all levels except for the last excited state and the J = 1 states of H^ and 

due to a higher numerical instability. In this paper, we extend our method to the 
J = 2 states and show that it remains advantageous with respect to the usual methods 
relying on Hylleraas coordinates. For convenience, we have regrouped all our results in 
the present paper ; all levels have been recomputed with the values of the most recent 
(2002) CODATA [12] ; and we also give their sensitivity with respect to the mass ratios. 
The accuracy for the J = 1 states of H^ and D^" has been improved to 10~ 14 a.u. On 
the whole, 57 vibration- rotation levels of Hj~, 79 levels of and 65 levels of HD + have 
been computed with a precision of metrological interest (10~ 12 a.u. or better). 

In the first section, we briefly recall the main features of our method (more details 
can be found in [1]). The second section is devoted to the presentation of numerical 
results. 



2. Method of resolution 



2.1. Hamiltonian 

Using centered Jacobi coordinates, the Hamiltonian of a three-body molecular ion with 
nuclear masses Mi and M 2 can be written as 



47re a V 2 2fi U \ A J 2// Q ||R/2 - r|| ||R/2 + r|| R 
where q is the electron charge, do the atomic Bohr radius, /X12 = MiM^/m^Mi + M 2 ) 
is the reduced mass of the two nuclei in units of m e , and l//xo = m e (l/M\ — 1/M 2 ). 
The dimensionless quantities R and r represent respectively the relative position of 
the two nuclei and the position of the electron with respect to their center of mass. 
The quantities P and p are the conjugate momenta. The term in p.P is the so-called 
symmetry breaking term, which vanishes for the homonuclear ions H^ and . 

2.2. Structure of the wave functions 

The rotational invariance yields the following separation between angular and radial 
variables: 

j 

v^(R,r) = D&rWA*) ®t M (R,P,0, (2) 
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where the radial coordinates R, p, ( have been defined in [I], if), 9, <f> are the Euler 
angles, and Df^ T are known angular functions related to the matrix elements of the 
rotation operator. As a result, for a given value of J and M, the wave functions are 
represented by 2 J + 1 unknown radial functions. 

Another symmetry of the Hamiltonian is the parity II. Since the parity only affects 
the angular part of the wave function, it is useful to introduce even and odd angular 
functions ; from the 2J+1 functions Dj% T one obtains J functions of parity II = (— 1) J+1 
and J + 1 functions of parity II = (— 1) J . 

In this paper, we consider bound states that are supported by the first Born- 
Oppenheimer curve (lsa g ) with a total angular momentum J < 2. These states have 
the symmetries S e , P°, D e and the corresponding separation between angular and radial 
variables is written, in the case of M = states: 

S e states: ^ 00 (R, r) = $°(R, p, () 

P° states: V W (R, r) = D x *(^, 9, 0) ®°(R, p, () 

H $ WftCJ 

= cos 9 <$>°(R,p,()-sm9 cos0 ^(R,p,C) 
D e states: * 20 (R, r) = Z$(V>, 9, (j)) $°(R, p, C) 

Dl*_ 1 (9,^,4>)-Dl* 1 (9,i>,<j>) 



V2 

Z^ 2 (^,0) + D o 2 *(£,?M) 



= 3 C ° S ^ ~ 1 * a (i2,P,C) - ^ sin(20) cos0 S 1 ^,, 

+ ^? siri 2 0cos(2 0) $ 2 (i2, p ,C) 
As a result, the radial part of the wave function is represented by J + 1 radial functions 

2.3. Factorization of the radial part 

The separation introduced in the previous section allows us to write down an effective 
Schrodinger equation for the radial part of the wave functions (for J > it is in fact a 
system of coupled equations involving the J + 1 radial functions $*(-R, p, £)). In order to 
regularize the divergence introduced by the centrifugal terms (for J > 0), it is necessary 
to factorize the radial wave functions. A general method of factorization was derived in 
One obtains, for P° states : 

®°(R,pX)= (c + ~) f(R,p,() G(R,p,C) 
& (R, p,() = P F(R, p,C) + P G(R, p, C) (6) 
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and for D e states : 



R\ 2 p 2 



+ 



C + 



^42 



F(R,p,C) + 
H{R,p,() 



c- 



R 



P_ 
2 



G(R,p, 



$ 1 ( J R,p,C) = VSpU + ^j F(R,pX) + V3p(c-^] G(R, 



P,0 



® 2 (R,P,0 



+ V3p(H(R,p,() 

V 3 F(R, p, C) + 4^ G(R, p, C) + 4^ P> ( 



2 v " ' " 2 2 
Using these expressions, we can write down a set of effective Schrodinger equations for 
the radial functions F,G,H appearing in ([MZ]), in which centrifugal terms of the type 
1/R 2 have disappeared. The effective Hamiltonian i? e // appearing in these equations 
can be found in [TJ for P° states. 



2.4- Exchange symmetry 

In the cases of and Dj, we have an additional symmetry corresponding to the 
exchange of the two nuclei. The wave functions are either symmetric or antisymmetric 
with respect to the exchange operator P\i- Like in the atomic case, we will note 
here spatially symmetric (respectively antisymmetric) as singlets (respectively triplets). 
The bound states of and D^" that are considered in this paper have the following 
symmetries : 1 S e , 3 P° and 1 D e . 

The effect of the exchange operator on radial functions is the transformation 
( — > —(. Thus the radial functions of singlet and triplet states differ by their behaviour, 
either symmetric or antisymmetric under the transformation ( — > — (. For the states 
considered here, the radial functions have the following properties : 

1 S e states: F = F (8) 
3 P° states: G = - F (9) 
x D e states: G = F and H = H (10) 

where F(R, p, C) = F(R, p, —(,). 



2.5. Numerical implementation 

Even though the centrifugal terms have been eliminated, divergences in l/r±, \jr% 
and 1/R remain due to the Coulomb potential {r\(2) = ||R ± r /2|| are the distances 
between the electron and the two nuclei) . These divergences can be regularized through 
multiplication of the Schrodinger equation by r 1 r 2 R. The Schrodinger equation is turned 
into a generalized eigenvalue problem, which is written as : 

A\^) = EB\^) (11) 
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with B = rir 2 R and A = BH e ff. 

Our numerical method to solve this problem has been explained in detail in [lj. It 
relies on the use of perimetric radial coordinates, defined by: 

x = ri + r 2 — R 

y = r 1 -r 2 + R (12) 
z — - 7*1 + r 2 + R, 
and Sturmian basis functions in the x, y, z coordinates 



n~ n 



P 4) = |<> ® |<> ® |nf>, (13) 



xi yi z I \ x I ^ \ y 

where |rt°) represents the function 

$ B (ou) = (u|n a ) = (-l) n v^4°) (aw) e- Qu/2 . (14) 

n is a non-negative integer and Ln the generalized Laguerre polynomials. a~ l is a 
length scale in the x direction, in the y and z directions. 

For and Dj, due to the exchange symmetry, some of the radial wave functions 
can be either symmetric or antisymmetric with respect to the transformation ( — ► — £, 
that is the exchange of y and z. In such cases we use a symmetrized or antisymmetrized 
basis : 

\i^xi 'i"y) ii'zl — V / 

To perform the numerical calculations, the basis is truncated at n x + n y + n z < N and 
n x < with ^ < N. If the basis is symmetrized (resp. antisymmetrized) we add 
the condition n y < n z (resp. n y < n z ) which reduces the size of the basis by a factor 
of about 2. Since the radial part of the wave function is represented by J + 1 radial 
functions, the size of the matrices representing the A and B operators also depends on 
J. The different cases are summarized in Table [U 

Because of their structure, all the terms in the Hamiltonian have strict coupling 
rules. A and B are then sparse band matrices having exactly the same shape ; the order 
of the basis vectors is chosen in order to minimize their width around the diagonal. 
The coupling rules and width are also reported in Table [1] the fact that the number of 
coupling rules increases with J is due to the factorization, which involves polynomials 
of degree J (see Eqs. ()6])-(J7j)). As a result, with increasing J the matrices become both 
larger and wider. 

The analytical calculation of the matrix elements of the various contributions to the 
Hamiltonian has been performed using the symbolic calculation language Mathematica 
4. The results are directly output in double precision FORTRAN code. The generalized 
eigenvalue problem is then diagonalized using the Lanczos algorithm. That gives the 
eigenvectors and eigenvalues in the energy range of interest. 
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J = 1 


J = 2 






HD+ 




HD+ 




HD+ 


radial wave function 


F = F 


F 


F 


F,G 


F, H = H 


F, G, 


basis size 


~ N tot /2 


Ntot 




2N tot 




3iV to t 


coupling rules 


|A7i Xj y j2 | — 2, 

|An x | + An, + An. | < 3 


An m | + |An v | + |An*| < 5 


lAna-^l < 6, 

| Ansel + |Any| + |An^| < 7 


number of coupling rules 


57 


450 


1707 


for N = 60 and JV^ = 15 


basis size 


11964 


23496 


23496 


46992 


35460 


70488 


width 


834 


1364 


3031 


4552 


6422 


8656 



Table 1. For all the computed levels, we have indicated in this table : the radial 
functions representing the wave function, and (if necessary) their symmetry with 
respect to the exchange of y and z ; the basis size as a function of N to t , where N to t is 
the number of Sturmian functions \n x , n y , n z ) verifying n x + n y +n z < N and n x < N x 
; the coupling rules between \n x , n y , n z ) and \n x + An x , n y + An y , n z + An z ) and the 
number of coupling rules. Finally, for a typical value of N and N x we give the basis 
size and the width of the matrix. 

3. Numerical results 

The energy levels of , and HD + are given in Tables [2l H] and [61 The mass 
ratios are taken from the 2002 CODATA [12]: M P /m e = 1836.15267261 and M D /m e = 
3670.4829652. The atomic unit of energy is 219474.6313705 cm -1 . All the digits shown 
in these tables are converged ; an accuracy of 10 -14 atomic unit, limited by the numerical 
noise, is achieved for most levels. Let us stress that these results were obtained with 
quite reasonable computation resources, i.e. a single standard workstation with 8 Go 
memory and double precision arithmetic. 

When the accuracy is sufficient to be sensitive to the mass ratios, we have computed 
the normalized sensitivity XdE/dX of the energy levels to the variation of the relevant 
mass ratio(s), that is the electron/nucleus mass ratio A = m e /Mp for or HD + 
(A = m e /Mo for Dj) and the ratio of the nuclear masses /x = Mp/Mp for HD + , 
following the method described in [3J. The sensitivities are given in Tables [31 [5] and [71 
Again, all the digits shown in these tables are converged ; they are in full agreement 
with the values published in [3] . From these tables it is possible to calculate the energy 
levels for any value of the mass ratios within their present area of uncertainty (in the 
2002 CODATA the relative uncertainty is 4.6 10" 10 for m e /M P} 4.8 10" 10 for m e /M D 
and 2 10~ 10 for Mp/Mp) while maintaining the same accuracy. 

The precision for the J = 1 states of and has been improved from 10~ n to 
10~ 14 a.u. with respect to earlier results pQ. As discussed in pQ, the increased numerical 
noise was due to round off errors that accumulate more rapidly during the Lanczos 
steps, because the matrices are ill-conditioned. However, for a given eigenvalue E n the 
accuracy of the result can be improved by solving the problem (A — XB) = E B \^>) 
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with A close to E n . In this way, the eigenvalue is converged during the first few Lanczos 
steps and the round off errors are greatly reduced. Of course, the downside is that a 
new matrix has to be diagonalized for each level, which multiplies the computation time 
by the number of computed levels. 

For J = 2 levels, the conditioning properties of the matrices are even worse than 
for J = 1 levels. It is not surprising that the conditioning should worsen when J 
increases, because their matrix elements are polynomials in n x , n y , n z , the degree of 
which increases with J since the degree of the polynomials appearing in the factorization 
increases (see expressions (Q, ([7])). However, for most levels it is still possible to achieve 
an accuracy of 10~ 14 a.u. Only for the last few levels (for example v — 16 and v — 17 
in H^, or v = 20 to 24 in D^") does the numerical accuracy drop to 10~ 13 , because of 
the increased basis size. On the whole, for most levels the precision has been improved 
by 4 or 5 orders of magnitude with respect to earlier works [U El E] except for the first 
vibrational levels of HD + (v = to v — 4) which had been computed with 10~ 14 or 
1CT 15 accuracy in [3]. We have checked that our results are in full agreement with those 
of that reference. 

Finally, let us note that for J = and J = 1, only the last vibrational level (or the 
last two levels) could not be converged to 10 -14 accuracy due to the limit in memory 
size. J = 2 levels demand more memory (since the matrices are both larger and wider, 
see Table [1]) which is why the level of convergence drops for the last three or four 
vibrational levels. Whenever the accuracy of our calculations did not reach 10 _9 -10 -10 
a.u. we have reported the more precise results obtained by R. E. Moss jUEJE]; although 
these results have been obtained with the 1986 CODATA values of the mass ratios, all 
the digits remain valid because the sensitivity of the last vibrational levels is very low. 

All these elements make it clear that our method becomes less and less convenient 
when J increases. It may be possible to extend the calculations to J = 3 states, but 
then two problems will occur: 

- the limit in memory will prevent the calculation of a larger number of levels ; 

- conditioning problems will become more severe ; unless an efficient method of pre- 
conditioning of the matrices can be devised, quadruple precision arithmetic is likely 
to be compulsory in this case which imposes new constraints in terms of memory and 
computation time. 

4. Conclusion 

We have shown that the use of perimetric coordinates allows to obtain high accuracy 
results for the nonrelativistic rotation- vibration energies of the hydrogen molecular ion 
and its isotopes, up to J = 2 ; our results suggest that this is probably the last value 
of J for which this method is advantageous with respect to usual methods relying on 
Hylleraas coordinates. The very accurate wavefunctions that we obtain enable the 
calculation of relativistic and radiative corrections with an improved precision, which is 
in progress. These results extend the range of optical transitions on which high-precision 
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measurements can be compared to theoretical predictions for metrological applications. 
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Table 2. Energies of the 1 S e , 3 P° and 1 D e bound levels of the molecular ion 
below the first dissociation limit, in atomic units. The star indicates that there is 
no bound level. a result taken from Ref. [5]. The first dissociation limit of is 
-1/2 (1 + me/Mp)- 1 = -0.499 727 839 712 26 a.u. 
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Table 3. Derivatives 10 2 XdE/dX of the energies of with respect to the 
electron/proton mass ratio A (in atomic units). The hyphens correspond to levels 
for which the accuracy of the calculation is not sufficient to be sensitive to a change of 
the recommended value of the proton/electron mass ratio. 



High accuracy results for the energy levels of the molecular ions H\ , B\ and HD + , up to J = 211 
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Table 4. Same as Table [21 for the D^" molecular ion. a result taken from Ref. [6]. 
The first dissociation limit of D+ is -0.499 863 815 247 21 a.u. 



High accuracy results for the energy levels of the molecular ions H^ > As an< ^ HD + , up to J = 212 
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Table 5. Derivatives 10 2 XdE/dX of the energies of with respect to the 
electron/deuteron mass ratio A (in atomic units). The hyphens correspond to levels 
for which the accuracy of the calculation is not sufficient to be sensitive to a change of 
the recommended value of the deuteron/electron mass ratio. 



High accuracy results for the energy levels of the molecular ions iij , B\ and HD + , up to J = 213 
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Table 6. Energies of the S e , P° and D e bound levels of the HD + molecular ion. 
The star indicates that there is no bound level, "result taken from Ref.[4]. The two 
dissociation limits of HD + are HD + — >D+H + (given in the caption of Table 0|) and 
HD+ — >H+D+ (given in the caption of Table [2]). 



High accuracy results for the energy levels of the molecular ions HJ > As an< ^ HD + , up to J = 214 
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Table 7. Derivatives of the energies of HD+ with respect to the electron/proton 
and proton/deuteron mass ratios, respectively noted A and u(in atomic units). The 
hyphens correspond to levels for which the accuracy of the calculation is not sufficient 
to be sensitive to a change of the recommended value of the mass ratios. 



